% Copyright (C) 2009 Indrek Mandre <indrek(at)mare.ee>

% usage: [S] = lsexp_min (D, A, X, COVM)
%
% Calculates the minimum squared residuals at exponent D.
%
% Arguments:
%   D - scalar, the exponent at which to find the minimum
%   A - N x 1 column matrix of base rows
%   X - N x M matrix of data rows
%   COVM - M x M x N covariance matrixes between data elements
%
% Returns:
%   S - minimized squared residuals for exponent D
%
% Required Octave-Forge packages: statistics
%

function [S] = lsexp_min (D, A, X, COVM)
  if abs(D) > 10
    S = 0;
    return;
  end
  N = size(X, 1);
  M = size(X, 2);
  C = ones(M, 1); % preliminary coefficients, a guess
  % we iterate up to 10 times to get the correct coefficients
  Y = A .^ D;
  for j=1:10
    mX = X;
    mY = Y;
    for i=1:N
      s = sqrt(sum(sum((C * C') .* COVM(:,:,i))));
      mY(i) = mY(i) ./ s;
      mX(i,:) = mX(i,:) ./ s;
    end
    oldc = C;
    %[C, sigma, r] = ols(mY,mX); % octave's ols() is 4x faster but gives us one less digit
    [C,bint,r] = regress (mY,mX);
    % iterate only as long as the difference in C-s is very small
    if (C - oldc)' * (C - oldc) < 1e-20
      break
    end
  end
  S = r' * r;
end

